function calcCellParamsFromNegVstep
% calcCellParamsFromNegVstep
%
% Calculate Rs, Ri, Cm, and tau from a trace with one or more 
% negative-going V steps
%
% NB: tailored for specific 4-step pulse with following parameters:
% sampleRate = 10000
% amplitude is extracted; typically -5
% squarePulseTrainNumber extracted, e.g. = 4
% squarePulseTrainISI, e.g. 0.2
% squarePulseTrainDelay, e.g. 0.1
% squarePulseTrainWidth, e.g. 0.1
% Also note that the last 20 msec are assumed to be usable as baseline!

% gs 2008 03 08
% ----------------------------------

% trace data
data = get(gcf, 'UserData');
currentTrace = data.currentTrace;
traw = currentTrace.data.ephys.trace_1;

% TODO: add some filtering options

% stimulus parameters
p = currentTrace.header.ephys.ephys.pulseParameters{:};
if isempty(p) % KLUGE!! 
    p = currentTrace.header.ephys.ephys.pulseParameters{2};
end
sr = p.sampleRate;
V = p.amplitude;
numsteps = p.squarePulseTrainNumber;
isi = p.squarePulseTrainISI;
delay = p.squarePulseTrainDelay;
dur = p.squarePulseTrainWidth;

% calculate baseline and baseline-subtract the whole trace
snippetAtStart = traw(1 : delay * sr - 1);
snippetAtEnd = traw((end - round(0.02 * sr)) : end); % last 20 msec
baseline = mean([snippetAtStart; snippetAtEnd]);
tbs = traw - baseline;

% calculate an average response
for n = 1 : numsteps
    startInd = round((delay + (n-1) * isi) * sr);
    endInd = startInd + round(dur * sr);
    tracelet(:,n) = tbs(startInd : endInd);
end
tmeanbs = mean(tracelet,2);

% make a corresponding time trace
x = 0:length(tmeanbs)-1;
x = x / sr;

% find peak, within first 10 msec
[Ipeak, peakX] = min(tmeanbs(1:round(0.01*sr)));
% assumes pulse is hyperpolarizing; response is downward

% steady state: average over last 20 msec of response
steadystate = mean(tmeanbs(end-round(0.02*sr):end));

% calculate Rseries and Rinput
Rs = 1000 * V / Ipeak;
Ri = 1000 * V / steadystate - Rs;

% steadystate-subtracted mean tracelet:
tmeanss = tmeanbs - steadystate;

% calculate Cm, based on time for signal to fall to 1/e of peak
efold = find(tmeanss(peakX:end) >= (Ipeak-steadystate)/exp(1));
dropInd = efold(1);
tau = x(dropInd);
Cm = 10^6 * (Rs + Ri) * tau / (Rs * Ri);

% graphics
figure('Color', 'w'); 
subplot(2,1,1)
plot(x,tmeanbs, 'b-')

subplot(2,1,2)
plot(x(peakX:end), tmeanss(peakX:end));
hold on
plot(tau, tmeanss(dropInd),'ro');

% command line display
disp(' ');
disp('Rs   Ri   Cm   tau');
disp([num2str(round(Rs)) '  ' num2str(round(Ri)) '  ' ...
    num2str(round(Cm)) '  ' num2str((tau))]);

% results -> workspace
cellParameters.Rs = Rs;
cellParameters.Ri = Ri;
cellParameters.Cm = Cm;
cellParameters.tau = tau;
assignin('base', 'cellParameters', cellParameters); 
